/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <pthread.h>

#include "estep.h"
#include "normalize.c"
#include "prating.c"

static const float log1mepsilon = -0.000010000050000333336;/* log (0.99999) */
static const float logepsilon = -11.512925464970228;       /* log (0.00001) */

static float
estep_thread (unsigned int               thread_no,
              unsigned int               n_threads,
              const float*               loga,
              unsigned int               rank,
              const float*               logbeta,
              unsigned int               n_images,
              const float*               c,
              unsigned int               n_labels,
              const float*               d,
              float*                     pz,
              const float*               priorz,
              const int*                 clamp,
              Rating*                    ratings,
              unsigned int               n_ratings)
{
  float logpriorz[n_labels];
  float log_n_labels_m1 = (n_labels < 2) ? 0 : log (n_labels - 1);

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      logpriorz[i] = log (priorz[i]);
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      if (i % n_threads == thread_no)
        {
          if (clamp[i] < 0 || n_labels < 2)
            {
              memcpy (pz + i * n_labels, logpriorz, n_labels * sizeof (pz[0]));
            }
          else
            {
              for (unsigned int j = 0; j < n_labels; ++j)
                {
                  if (((unsigned int) clamp[i]) == j)
                    {
                      pz[i * n_labels + j] = log1mepsilon;
                    }
                  else
                    {
                      pz[i * n_labels + j] = logepsilon - log_n_labels_m1;
                    }
                }
            }
        }
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          if (ratings[i].image % n_threads == thread_no)
            {
#define logafunc(n) loga[ratings[i].worker * rank + n]
#define cfunc(n) c[n_labels * n + j]
#define dfunc(n, k) ((k == j) ? 0 : d[n_labels * n + k])

              pz[j + ratings[i].image * n_labels] += 
                logprating (ratings[i].rating,
                            logafunc,
                            rank,
                            logbeta[ratings[i].image],
                            cfunc,
                            dfunc,
                            n_labels);
#undef dfunc
#undef cfunc
#undef logafunc
            }
        }
    }

  float q = 0;

  for (unsigned int i = 0; i < n_images; ++i)
    {
      if (i % n_threads == thread_no)
        {
          q += safe_exponentiate_normalize (pz + i * n_labels, n_labels);
        }
    }

  return q;
}

typedef struct _EstepData EstepData;
struct _EstepData 
{
  unsigned int               n_threads;
  const float*               loga;
  unsigned int               rank;
  const float*               logbeta;
  unsigned int               n_images;
  const float*               c;
  unsigned int               n_labels;
  const float*               d;
  float*                     pz;
  const float*               priorz;
  const int*                 clamp;
  Rating*                    ratings;
  unsigned int               n_ratings;
};

typedef struct _EstepArg EstepArg;
struct _EstepArg 
{
  EstepData* data;
  unsigned int thread_no;
  float q;
};

static void*
estep_start (void* arg)
{
  EstepArg* estep_arg = (EstepArg*) arg;

  estep_arg->q = estep_thread (estep_arg->thread_no,
                               estep_arg->data->n_threads, 
                               estep_arg->data->loga, 
                               estep_arg->data->rank, 
                               estep_arg->data->logbeta,
                               estep_arg->data->n_images, 
                               estep_arg->data->c, 
                               estep_arg->data->n_labels,
                               estep_arg->data->d,
                               estep_arg->data->pz,
                               estep_arg->data->priorz,
                               estep_arg->data->clamp,
                               estep_arg->data->ratings,
                               estep_arg->data->n_ratings);

  return NULL;
}

float
estep (unsigned int               n_threads,
       const float*               loga,
       unsigned int               rank,
       const float*               logbeta,
       unsigned int               n_images,
       const float*               c,
       unsigned int               n_labels,
       const float*               d,
       float*                     pz,
       const float*               priorz,
       const int*                 clamp,
       Rating*                    ratings,
       unsigned int               n_ratings)
{
  pthread_t threads[n_threads];
  EstepArg args[n_threads];
  EstepData data = { .n_threads = n_threads,
                     .loga = loga,
                     .rank = rank,
                     .logbeta = logbeta,
                     .n_images = n_images,
                     .c = c,
                     .n_labels = n_labels,
                     .d = d,
                     .pz = pz,
                     .priorz = priorz,
                     .clamp = clamp,
                     .ratings = ratings,
                     .n_ratings = n_ratings };

  for (unsigned int k = 1; k < n_threads; ++k)
    {
      args[k].data = &data;
      args[k].thread_no = k;
      pthread_create (threads + k, NULL, estep_start, args + k);
    }

  float q = estep_thread (0, n_threads, loga, rank, logbeta,
                          n_images, c, n_labels, d, pz, priorz,
                          clamp, ratings, n_ratings);

  for (unsigned int k = 1; k < n_threads; ++k)
    {
      pthread_join (threads[k], NULL);
      q += args[k].q;
    }

  return q;
}
